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1  Introduction 

The  nonlinear  transverse  dynamics  of  single  longitudinal  mode  lasers  in  cavities  is 
described  by  the  Maxwell-Bloch  (MB)  equations  that  account  for  the  interaction  of 
the  electromagnetic  field  with  the  active  media  inside  the  cavity.  Several  different 
order  parameter  equations  have  been  derived  from  the  MB  equations  and  are  cur¬ 
rently  used  to  gain  some  understanding  of  transverse  pattern  forming  instabilities  in 
lasers  (see,  e.g.,  the  recent  reviews  [1]  and  [2]).  The  study  of  these  order  parameter 
equations  (OP)  is  interesting  not  only  for  the  analysis  of  the  transverse  nonlinear 
dynamics  of  lasers,  but  also  from  the  more  general  point  of  view  of  pattern  forma¬ 
tion  in  general  spatially  extended  systems,  because  the  same  OP  equations  can  be 
used  to  describe  other  physical  systems  that  exhibit  the  same  basic  destabilization 
mechanism  [3]. 

The  derivation  of  the  OP  equations  from  the  MB  equations  is  based  on  the  slow 
envelope  assumption,  i.e.,  on  the  assumption  that  the  amplitudes  of  the  physical 
fields  are  small  and  depend  slowly  on  time  and  on  the  transversal  spatial  scales. 
But  the  OP  equations  are  commonly  applied  without  observing  the  slow  envelope 
requirement,  and  it  is  not  infrequent  to  find  in  the  literature  solutions  of  the  OP 
equations  that  do  not  exhibit  only  long  scales.  This  is  normally  justified  by  saying 
that  the  OP  equations  retain  the  essential  dynamics  of  the  system  and  thus  its  range 
of  application  can  be  somehow  extended  and  they  still  provide  ’’qualitative”  infor¬ 
mation  about  the  patterns  and  instabilities  that  take  place  in  the  original  physical 
system  (see  [1,  2]  and  the  references  therein). 

The  main  goal  of  this  research  work  is  to  go  beyond  the  qualitative  description 
and  quantitatively  validate  the  OP  equations  to  show  how  accurate  and  precise  are 
the  predictions  they  provide  about  the  dynamic  states  and  transitions  that  take 
place  in  the  transverse  laser  profiles. 

To  this  end,  we  will  first  apply  multiple  scales  techniques  to  the  MB  equations  to 
derive  the  corresponding  OP  equation  and  its  coefficients  in  terms  of  the  physical 
parameters  of  the  problem.  For  the  laser  parameter  values  that  we  will  consider, 
namely  class  C  with  small  detuning,  the  resulting  OP  equation  is  a  complex  Swift- 
Hohenberg  (CSH)  equation  [4,  5].  In  this  case,  due  to  the  small  value  of  the  detuning, 
the  diffusion  of  the  system  acts  on  a  much  slower  time  scale  than  the  dispersion 
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and,  if  the  slow  envelope  assumption  is  consistently  applied,  this  scale  disparity  is 
inherited  by  the  CSH  equation  that  necessarily  includes  terms  of  different  asymptotic 
order  and,  in  in  some  particular  regimes,  it  even  allow  us  to  further  simplify  the 
resulting  CSH  to  a  complex  Ginzburg-Landau  (CGL)  equation.  This  unavoidable 
asymptotic  nonuniformity  of  the  resulting  CSH  equation  is  not  properly  taken  into 
account  in  the  usual  analyses  in  the  literature,  where,  invoking  the  argument  of 
the  qualitative  only  scope  of  the  OP  description,  a  CSH  equation  with  all  terms  of 
the  same  order  is  normally  used  to  analyze  this  laser  instability.  In  order  to  clarify 
this  point,  we  carry  out  in  this  research  effort  numerical  integrations  of  both,  the 
MB  equations  in  an  extended  spatial  domain  and  the  CSH  equation,  that  allow  us 
to  confirm  that  the  quantitatively  correct  weakly  nonlinear  description  of  the  MB 
system  is  given  by  the  above  mentioned,  asymptotically  nonuniform  CSH  equation 
(and,  in  some  cases,  by  the  simpler  CGL  equation). 

This  report  is  organized  as  follows.  In  the  next  Section  we  introduce  our  physical 
model  for  the  laser  transverse  dynamics,  the  MB  equations,  and  study  the  basic 
linear  stability  characteristics  of  the  non-lasing  state  in  a  class  C  laser  with  small 
detuning.  The  derivation  of  the  corresponding  CSH  equation  that  describe  the  near 
onset  dynamics  of  the  system  is  carried  out  in  Section  3  and  the  simplest  pat¬ 
terns  (travelling  waves  and  spiral  waves)  displayed  by  the  system  and  their  stability 
properties  are  analyzed  in  Section  4.  Finally,  Section  5  contains  some  numerical 
simulations  of  the  CSH  equation  performed  in  order  to  check  the  presence  of  two 
asymptotically  different  scales,  and  the  description  and  results  of  the  procedure  im¬ 
plemented  in  order  to  quantitative  compare  the  solutions  the  MB  and  CSH,  together 
with  a  brief  explanation  of  the  numerical  methods  used  to  integrate  both,  the  MB 
and  CSH  equations. 
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2  Maxwell-Bloch  model 


The  dynamics  of  a  two-level  laser  in  the  2D  plane  section  transverse  to  the  main 
propagation  direction  of  the  electromagnetic  wave  is  governed  by  the  well-known 
Maxwell-Bloch  (MB)  equations  [6]  that,  in  the  simple  case  of  a  single  longitudinal 
mode,  can  be  written  as 


dE 

~dt 
dP  _ 

~dt  ~ 
d  N 

~dt  ' 


-  ia'V2E  —  crE  +  aP, 

(1  +  ito)P  +  (r  -  N)E, 
-bN+^(EP  +  EP). 


(1) 

(2) 

(3) 


We  are  using  here  the  same  nondimensional  formulation  as  in  [7],  where  E(x,y,t) 
and  P(x,y,t )  represent  the  complex  electric  and  polarization  fields,  N(x,y,t )  is  the 
real  valued  field  of  the  population  inversion,  a  >  0  measures  the  strength  of  the 
diffraction,  a  >  0  and  Q  account  for  the  cavity  losses  and  the  cavity  detuning, 
r  is  the  pumping  parameter,  b  >  0  is  the  decay  rate  of  the  population  inversion, 
V2  =  d2/dx 2  +  d2/dy2  is  the  2D  laplacian  operator  and  the  bar  stands  for  the 
complex  conjugate.  In  the  equations  above  the  size  of  the  fields  have  been  scaled  in 
order  to  remove  the  coefficients  of  the  nonlinear  terms  and  time  is  measured  using 
the  characteristic  decay  time  of  the  polarization.  Note  also  that  the  diffraction 
coefficient,  a,  can  be  set  to  1  without  loss  of  generality  (as  we  will  do  hereafter)  just 
by  scaling  the  space  variables  x  =  (x,  y )  with  the  diffraction  length  y/a. 

Lasers  are  usually  classified  according  to  the  relative  value  of  the  decay  rate  of 
the  population  difference  yy,  the  cavity  decay  rate  n,  and  the  decay  rate  of  the 
polarization  7j_:  for  a  class  A  laser,  7j_,  yy  S>  K]  for  a  class  B  laser,  yj_  3>  k  yy ; 
and  all  decay  rates  are  comparable  for  a  class  C  laser,  n  ~  yj_  ~  yy,  see  e.g.  [8].  We 
will  focus  our  attention  on  class  C  lasers  that,  in  the  nondimensional  formulation 
of  eqs.  (l)-(3)  where  time  has  been  scaled  with  y_i_  and  a  =  ft/y_i_  and  b  =  yy/y_L, 
correspond  to  a  ~  1  and  b  ~  1. 

The  linear  stability  properties  of  the  basic  nonlasing  state  ( E ,  P ,  N )  =  (0,  0,  0) 


3 


Figure  1:  MB  zero  state  neutral  stability  curves  for  Cl  <  0,  Cl  =  0  and  Cl  >  0  (shading  indicates 
instability). 

are  readily  calculated  inserting  the  ansatz 

E(x,y,t) 

P(x,y,t) 

_  N(x,y,t) 

into  the  linearized  version  of  the  MB  and  solving  the  resulting  eigenvalue  problem, 


A  +  cr  +  i  k2  —a  0 

'  E  ' 

"  0  ' 

— r  A  +  1  +  iO  0 

P 

= 

0 

0  0  A  +  b  _ 

N 

k 

_  0  _ 

which  gives  the  following  stability  exponents,  A,  as  a  function  of  the  wavevector, 
k  ( kx ,  ky ) , 


E 

P 

N 


_ik-x+A£ 


J  lr 


A  =  -b,  (4) 

cr  +  1  +  i (k2  -  Q)  \/((rr  -  1)  —  i {k2  -  fl))2  +  4rcr 
A  -  g  ±  2  ’  (5) 

with  k2  =  |  k  | J .  After  some  tedious  algebra,  one  can  see  that  the  zero  solution  is 
stable  if  the  pump,  r,  is  below  the  critical  value 


rc 


, ,  (k2-ny 

(a  +  l)2  ’ 


and  unstable  otherwise  [4,  5]. 
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Re(A) 


n  <  o 


Re(A) 


n  >  o 


Figure  2:  Growth  rate  of  the  MB  almost  neutral  modes  near  onset  for  O  <  0  and  fl  >  0. 

Figure  1  shows  the  resulting  neutral  stability  curves,  i.e.,  the  critical  value  of  the 
pump  rc,  against  the  wavenumber  k  for  different  values  of  Q.  When  D  <  0,  the 
most  unstable  mode  is  k  —  0,  while  when  hi  >  0  the  most  unstable  modes  lie  along 
a  circumference  (in  the  plane  kx-ky)  of  radius  k  =  y/H.  These  are  the  modes  that 
first  become  unstable  when  the  pump  is  increased  over  the  lasing  threshold  rc.„Mn, 
which  is  given  by  rCimin  =  1  +  (g+~)2  for  hi  <  0  and  rC)mm  =  1  for  hi  >  0. 

Let  us  consider  now  that  the  system  is  close  to  threshold:  r  =  rcmm  +  e,  with 
e|  -C  1.  The  order  of  magnitude  of  the  unstable  wavenumbers  can  be  determined 
from  Fig.  2.  If  0  <  0,  then  k  ~  y/\t\  and,  using  scale  separation  and  applying 
solvability  conditions,  a  well  known  complex  Ginzburg-Landau  (CGL)  equation  is 
straightforwardly  obtained  [7]  for  the  weakly  nonlinear  description  of  the  system 
dynamics.  If,  on  the  other  hand,  O  >  0,  then  the  circumference  of  neutral  modes  in 
the  plane  kx-ky  becomes  a  thin  annulus  of  unstable  modes  of  width  A k  V\~e\  (see 
Fig.  2).  In  this  case  a  a  pair  of  counter-propagating  coupled  CGL  equations  can  be 
obtained  for  ID  and  quasi-ID  configurations  [7,  9,  10],  but  the  OP  equation  that 
completely  describes  this  incipient  destabilization  in  2D  is  not  currently  known,  and 
it  still  constitutes  one  of  the  open  problems  in  pattern  formation  [3]. 

This  report  covers  the  transitional  case  of  small  detuning  1 0|  <C  1  in  a  class  C 
laser,  which  is  described  by  a  complex  Swift-Hohenberg  (CSH)  equation  [4,  5].  The 
CSH  can  be  also  obtained  for  class  A  lasers,  and  even  for  class  B  lasers  (e.g.,  CO2 
and  semiconductor  lasers),  but  starting  from  a  more  specific  set  of  equations  [11,  12], 
A  detailed  derivation  of  the  CSH  from  the  MB  equations,  where  special  emphasis 
is  placed  on  the  asymptotic  order  of  its  different  terms,  is  presented  in  the  next 
section. 
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3  Order  parameter  equation  for  small  detuning 


In  this  section  we  derive  the  slow  envelope  equation  that  describes  the  transverse 
laser  patterns  displayed  by  the  MB  equations  (l)-(3)  right  after  onset,  in  the  small 
detuning  case. 

Near  the  destabilization  threshold,  the  weakly  nonlinear  dynamics  of  the  system 
is  dominated  by  its  linear  part  and  the  order  of  magnitude  of  the  slow  scales  that 
the  system  exhibits  can  be  obtained  from  the  study  of  the  behavior  of  the  branch 
of  eigenvalues  responsible  for  the  stability  change.  The  small  nonlinearities  of  the 
problem  do  not  play  any  role  in  the  slow  scale  selection;  they  just  define  the  charac¬ 
teristic  size  of  the  resulting  pattern  through  the  saturation  of  the  instability  growth, 
and  they  will  be  taken  into  account  later  on. 

For  small  detuning,  the  modes  of  the  system  that  first  become  unstable  are  those 
corresponding  to  small  wavevectors,  and  their  linear  stability  characteristics  are 
given  by  expression  (5)  with  the  +  sign,  which  can  be  expanded  as 


Re(A) 
Im  (A) 


a 


cr  +  1 


(r  -  1)  - 


cr 


(a  +  1)' 


:(k2 


k 2  + 

cr  +  1 


+  ... 


5 


^)2  H - j 


(6) 

(7) 


in  the  limit  of  small  detuning,  small  displacement  from  threshold  and  small  wavevec- 
tor,  i.e., 

\£l\  <C  1,  | r  —  1|  -C  1  and  fc«l.  (8) 


From  the  expression  of  the  growth  rate  (6)  we  can  see  that  the  characteristic  slow 
time  is  given  byr~l/|r  —  1|S>1  and  that  the  range  of  almost  neutral  modes  goes 
up  to  k  rsJ  |  V  -II1/4  (see  Fig.  3).  The  growth  rate  of  the  modes  with  !;>  |r  —  l]1/4  is 
much  larger  and  negative  and  therefore  they  are  quickly  damped  out  in  the  slow  time 
scale.  Note  that  in  the  scaling  arguments  above  we  are  implicitly  assuming  that  |0 1 
is,  at  most,  of  the  order  of  | r  —  1 1 4/2.  We  are  not  considering  the  case  1 0 1  |r  —  1 1 4/2 
and  Q  >  0  (Q  <  0  is  stable  and  far  from  threshold),  in  which  a  narrow  band  of 
almost  neutral  modes  of  size  A k  ~  |r  —  l)1/2/ -y/jTTf  develops  around  the  critical 
wavevector  circumference  of  radius  kc  ~  A k,  because  it  corresponds  to  a 

completely  different  destabilization  type,  namely  to  that  represented  in  Fig.  2  for 

Q  >  0. 
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r  —  ll1/4 


Figure  3:  Sketch  of  the  growth  rate  of  the  most  unstable  MB  modes  right  after  onset,  for  |0|  <C  1, 
\r  —  1|  <C  1  and  k  <C  1. 

In  other  words,  the  behavior  of  the  growth  rate  ensures  that,  for  large  times,  the 
system  can  only  exhibit  spatial  scales  that  are  of  the  order  1/| r  —  1| 1  /4  3>  1  or  larger, 
but  no  smaller  ones,  because  those  are  damped  out  by  the  dominant  linear  effects. 

On  the  other  hand,  the  associated  mode  frequencies  given  by  eq.  (7)  are  always 
large  as  compared  with  the  damping  terms  of  the  growth  rate.  This  indicates  that 
the  effect  of  dispersion  on  the  neutral  modes  of  the  system  is  faster  and  much  stronger 
than  that  of  diffusion.  The  envelope  equation  that  describes  the  slow  dynamics  of 
the  system  has  to  take  into  account  these  two  effects  simultaneously  and  therefore  it 
is  necessarily  going  to  contain  terms  of  different  asymptotic  orders,  which  can  never 
be  rescaled  to  be  of  the  same  order.  This  is  not  the  case  in  the  transition  represented 
in  Fig.  2  for  f!  <  0,  where  dispersion  and  diffusion  are  both  of  the  same  order  ~  k 2 
and  the  weakly  nonlinear  dynamics  of  the  system  is  described  by  an  asymptotically 
uniform  complex  Ginzburg-Landau  equation  [7]. 

The  need  to  retain  terms  of  different  orders  made  the  derivation  of  the  envelope 
equation  presented  in  [4,  5]  quite  involved.  They  first  assumed  a  particular  scaling  for 
the  fields  and  the  slow  spatial  and  temporal  variables,  then  they  proceeded  to  derive 
the  envelope  equation  and  its  next  order  correction,  and,  finally,  they  collapsed  back 
the  two  term  expansion  into  the  original  (unsealed)  variables  in  order  to  obtain  a 
single  envelope  equation. 

The  derivation  of  the  envelope  equation  that  we  present  below  is,  we  believe, 
much  more  simple  because  we  do  not  assume  a  priori  any  scaling  relation  between 
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the  variables,  we  just  look  for  small  amplitude  solutions  that  depend  slowly  on  space 
and  time  (similar  scaling  free  derivations  can  be  found  in  [13]  and  [14]  in  the  context 
of  convection  and  water  waves). 

More  precisely,  the  neutral  mode  at  threshold,  r  —  1  and  0  =  0,  corresponds  to 
A  =  0  and  k  —  0,  with  associated  eigenvector 


and  suggests  an  expansion  for  the  solution  of  the  MB  equations  (l)-(3)  of  the  form 

E(x,y,t) 

P(x,  y,  t )  =  V0i>  +  Vi  VV  H - ,  (9) 

.  N(x,y,t)  _ 

where  the  complex  amplitude  ij)(x,y,t)  has  to  verify  a  solvability  condition  that  is 
also  expanded  as 

ipt  =  aoV’  +  oqVV  H —  ,  (10) 

and  is  precisely  the  envelope  equation  that  we  are  looking  for.  The  expansions  above 
include  all  possible  combinations  of  the  powers  of  the  small  quantities 

\n\  «  1,  |r  —  1|  -C  1  and  |^|  <C  1,  |  V2^|  <  1,  |  VV|  -C  1, . . . 

that  are  treated  as  independent  parameters,  only  constrained  by  the  weakly  nonlin¬ 
ear  level  of  this  approach,  which  requires  essentially  that 

•  •  •  -C  |  V‘V|  -C  |  v2^|  -C  It/’l  1  and  •  •  •  -C  \^t\  -C  |^|  -C  1 . . . , 

i.e.,  small  envelope  amplitude  and  slow  time  and  space  dependence  (we  have  also 
anticipated  the  fact  that  only  even  order  spatial  derivatives  will  be  produced  because 
of  the  particular  structure  of  the  MB  equations). 

When  the  above  expansions  (9)  and  (10)  are  inserted  in  the  MB  equations  (1)- 
(3),  a  linear  nonhomogeneous  singular  system  is  obtained  at  each  order.  And,  as 
it  is  the  standard  procedure  in  multiple  scale  methods  to  ensure  that  no  secular 
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terms  are  present  [15,  16],  the  nonhomogeneous  terms  of  these  systems  must  satisfy 
a  solvability  condition  that  provides  the  contribution  of  this  order  to  the  envelope 
equation. 

The  linear  terms  of  the  envelope  equation  (10)  can  be  easily  anticipated  because 
they  correspond  to  the  Taylor  expansion  of  the  critical  eigenvalue  branch  A (k2,  r,  11) 
(eq.  (5)  with  the  +  sign)  at  (k2,r,  fl)  =  (0,1,0)  (see,  e.g.,  [3]), 


A(o,i,o)V’  + 


dX 

dr 


(r  -  1)^  + 


(0,1,0) 


<9A 


Qi/j  — 


d\ 


1  d2  X 


2  d(k 


2)2 


vV- 


d2  X 


(0,1,0) 


dfld{k2 


(0,1,0) 

2„ 


d(k 


V2?/>  + 


(0,1,0) 


+  -  — 
^  2  dvt2 


(0,1,0) 
Q2^  + 

(0,1,0) 


and  thus,  according  to  expressions  (6)  and  (7),  are  of  the  form 


a 

a  +  1 


((r 


1)  —  ill)'?/'  +  i 


(>  +  l) 


V2^ 


a 

( cr  +  l)3 


(fi  +  V2)2^  +  •  •  • 


(11) 


The  computation  of  the  nonlinear  terms  requires  to  take  the  expansions  (9)  and 
(10)  to  the  MB  equations,  which  should  be  first  rewritten  in  the  more  convenient 
form 


a 

—a 

0  " 

"  E  ' 

d 

~  ~dt 

’  E  " 

i  W2E 

0 

-1 

1 

0 

P 

P 

+ 

(r  -  1  )E  -  i VtP 

+ 

—NE 

0 

0 

b  _ 

N 

N 

0 

J  (EP  +  EP)_ 

(12) 


with  the  dominant  terms  grouped  in  the  left  hand  side.  The  first  order  nonlinear 
contribution  to  the  expansions  (9)  and  (10)  is  given  by 


Vl'i/j]2  and  a|^|2, 


where  the  vector  V  and  the  scalar  a  are  obtained  from 


a 

—cr 

0  " 

"  0  " 

-1 

1 

0 

V=  -aV0  + 

0 

0 

0 

b  _ 

1 

This  system  can  only  be  solved  if  the  right  hand  side  is  orthogonal  to  the  solution 
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of  the  adjoint  problem 


Va 

vo 


and  this  solvability  condition  yields 


1 

a 

0 


1/  = 


0 

0 

i 

b 


and  a  =  0. 


Thus  there  is  no  contribution  of  this  term  to  the  envelope  equation,  in  other  words, 
it  is  not  a  resonant  term.  At  next  order  we  have  the  terms 


Wp\p>\2  and  pp>\p\2, 


which  must  verify 


a 

— 0 

0  " 

0 

-1 

1 

0 

W=  -pv0  + 

1 

b 

0 

0 

b  _ 

0 

and,  again,  we  apply  the  solvability  condition  to  obtain  the  following  coefficient  of 
the  envelope  equation 

^  6(l  +  cr)‘ 


Finally,  after  collecting  the  leading  nonlinear  contribution  above  and  the  linear 
terms  in  (11),  the  envelope  equation  reads 


A  = 


a((r  —  1)  —  ifl) 
0+1 


ip  + 


’0+1) 


vV- 


a 


(u+l): 


:(fi- 


-V2)V 


0- 


6(o+l) 


V#|2  +  • '  •  ,  (13) 


which  is  the  well  known  complex  Swift-Hohenberg  equation  (CSH)  (see  [4,  5]).  It  is 
important  to  stress  that  this  equation  has  to  be  solved  in  the  limit 


\n\  -C  1  and  |r  —  1|  -C  1, 


(14) 


and  that  its  solutions  must  verify  the  slow  envelope  assumption 

•••  <  |v4^|  <  |v2^|  <  m  <  1. 


(15) 
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And  therefore,  as  it  was  expected  from  the  study  of  the  growth  rate,  the  CSH  above 
contains  different  order  terms:  the  effect  of  dispersion  is  always  much  bigger  than 
the  dissipation  of  the  double  diffusive  term. 

We  now  introduce  the  small  parameter  0<£<1  and  the  scaling  of  the  variables 


(r  -  1) 


G2 


(cr  +  l)2 


(a+1)2-cs\ 


cr* 


a 


(. x,y )—  ^  ^-^\x,y)£. |  ip  =  e  \/b— —  ^  <pe  ,Q  =  —  +  —  cue,  (16) 

v  <y  '  a  <7 


which  brings  eq.  (13),  after  dropping  tildes,  to  the  more  convenient  form 


<h  =  a<i>  +  iV2</>  -  (p\(p\2  -  2£tuV2(j)  -  £2V4</>,  (17) 


where  the  rescaled  pump,  a,  and  detuning,  u>,  are  now  order  one  parameters.  Note 
that  a  can  be  also  absorbed  by  means  of  a  trivial  change  of  variables,  but  we  prefer 
to  keep  it  because  later  on,  when  we  perform  some  numerical  integrations  of  eq. 
(17),  it  will  allow  us  to  make  the  scaled  size  of  the  domain  independent  of  e. 

In  the  scaled  CSH  eq.  (17)  the  dominant,  dispersion  induced  dynamics  is  of 
order  one:  the  typical  length  scale  exhibited  by  the  dispersive  patterns  (that  was 
~  1/ y/\r  —  1|  1  in  the  original  variables)  is  now  hdisp  ~  1.  On  the  other  hand, 

the  small  coefficient  multiplying  the  double  diffusive  term  indicates  that  the  system 
is  capable  of  displaying  scales  as  small  as  hdiff  ~  a fi  (  ~  1/ y/\r  —  1|  S>  1  in  the 
original  variables)  before  diffusion  becomes  overwhelming  and  everything  is  damped 
out,  recall  Fig.  3. 

The  most  simple  solutions  of  eq.  (17)  are  those  that  only  exhibit  scales  ~  <5disP  ~  1 
as  £  — >  0.  For  those  dispersive  states,  the  spatial  derivatives  remain  of  order  one  and 
thus  the  two  last  terms  in  eq.  (17)  represent,  in  the  limit  £  — >  0,  a  small  correction 
that  can  be  neglected  to  give 


4>t  =  a(f>-\-  iV20  -  (p\(j)\2.  (18) 

This  is  a  very  interesting  complex  Ginzburg-Landau  equation  (CGL)  [3,  17,  18]: 
purely  dispersive  but  with  a  dissipation  that  comes  from  the  nonlinear  term  only. 
This  CGL  appears  in  a  straightforward  way  from  the  small  detuning  MB  equations 
due  to  the  dominant  character  of  the  dispersion  but,  so  far,  we  have  not  seen  it 
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analyzed  in  the  literature.  Instead,  the  equations  that  are  currently  used  to  describe 
the  weakly  nonlinear  dynamics  of  the  MB  equations  in  the  small  detuning  case  are 
either  (i)  the  unsealed  CSH  eq.  (13),  but  integrated  without  observing  conditions 
(14)  and  (15),  see  [7,  4,  5,  19],  or  (ii)  a  CSH  eq.  scaled  similarly  to  eq.  (17),  see  [20], 
but  again  studied  for  e  ~  1,  when  £  -C  1  is  the  only  physically  relevant  regime  from 
the  point  of  view  of  the  slow  envelope  description  of  the  transverse  laser  pattern 
dynamics. 

Besides  lasers,  the  CSH  equation  has  been  used  as  a  model  for  other  non-linear 
optical  systems  such  as  optical  parametric  oscillator  [21,  22]  and  photorefractive 
materials,  where  the  equation  qualitatively  reproduces  experimental  results  [23]. 


12 


4  Patterns 


This  section  is  dedicated  to  the  study  of  some  basic  solutions  of  the  CGL  eq.  (18), 
namely:  spatially  uniform  travelling  waves  and  spirals.  Note  that  all  solutions  of 
the  CGL  are  approximate  solutions  (up  to  O(e)  corrections)  of  the  CSH,  and  once 
a  stable  solution  of  the  CGL  is  found,  in  order  to  be  sure  that  it  is  a  stable  state 
of  the  CSH,  one  must  also  check  its  stability  under  perturbations  containing  small 
diffusive  scales,  whose  dynamics  is  not  contained  in  the  CGL. 

As  a  previous  step,  we  study  the  stability  properties  of  the  nonlasing  solution 
(j)  =  0.  The  solutions  of  the  linearized  version  of  the  CSH  in  a  square  periodic 
domain  can  be  written  in  the  form 


</> 


r\j 


ikx+A  t 


where  A  is  given  by 

A  =  a  —  ik2  +  2  eujk2  —  £2kA , 

with  k  =  |  k  |  and  0  <  £  C  1.  Two  different  kinds  of  perturbations  can  be  distin¬ 
guished.  Those  with  wave  vector  k  ~  1,  i.e.,  those  exhibiting  only  dispersive  scales 
5disp  ~  1,  for  which  the  dispersion  relation  above  simplifies  to 

A  =  a  —  ik2  +  . . . ,  (19) 

and  thus  are  stable  for  a  <  0  and  become  unstable  for  a  >  0.  The  evolution  of 
these  perturbations  is  well  represented  by  the  CGL,  which  could  have  been  used  to 
obtain  their  stability  properties.  The  second  kind  of  perturbations  are  those  with 
very  high  wave  vector  k  ~  K/ y/e  3>  1  and  K  ~  1,  which  have  a  typical  wavelength 
of  the  order  of  the  small  diffusive  scales  Aiiff  ~  a J~£  -C  1,  and  whose  growth  rate  is 
given  by 

Re(A)  =  a  +  2 uK2  -  A'4.  (20) 

For  negative  detuning  uj  <  0  these  small  scale  perturbations  become  unstable  for 
positive  a,  but  for  u  >  0  the  destabilization  takes  place  for  a  >  —  to2  at  a  critical 
wavenumber  k  =  y/u/ \fe  1.  The  dynamics  of  these  perturbations  is  not  included 
in  the  CGL  and  its  analysis  requires  to  deal  with  the  complete  CSH. 

The  structure  of  the  growth  rate  of  these  perturbations  is  sketched  in  Fig.  4  and, 
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Re(A) 


CGL 


CSH 


Figure  4:  Growth  rate  of  the  perturbations  of  the  zero  solution  of  the  CSH  eq.  (17)  for  u>  >  0  and 
to  <  0,  with  the  two  wavenumber  regions  shown:  k  ~  1  (dispersive  perturbations)  and  k  ~  1  /y/e 
1  (diffusive  perturbations). 

as  it  was  expected,  it  is  just  a  rescaled  version  of  that  obtained  in  the  previous  section 
for  the  MB,  see  Fig.  3.  The  resulting  stability  region  of  the  nonlasing  solution  is  the 
dark  shaded  area  in  the  a  —  u  plot  in  Fig.  5. 

Inside  the  linear  stability  region  depicted  in  Fig.  5,  the  zero  solution  is,  in  fact, 
globally  stable,  that  is,  all  solutions  of  the  CSH  for  all  possible  initial  conditions 
decay  to  0  as  time  increases.  In  order  to  see  this,  we  multiply  eq.  (17)  by  0,  add 
the  complex  conjugate  and  integrate  over  the  domain  D  to  obtain 

4  f  I0|2  =  2  /  !0|2(«  -  10!') +  i  f  (0V20-0V20) 

Gu  J  D  J  D  J  D 

-2au  [  (0V20  +  0V20)  +  e2  f  (0V40  +  0V40). 

J  D  JD 

If  we  now  take  into  account  that 

0V20  =  V  •  (0V0)  -  |V0|2, 

and  apply  Green’s  formula,  the  second  integral  above  can  be  expressed  as 

=  o. 

that  is  zero  because  of  the  periodicity  boundary  conditions  at  the  boundary  5D  of 
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Figure  5:  Stability  diagram  for  the  zero  solution  and  the  TW.  Light  shading:  TW  stable.  Dark 
shading:  zero  solution  stable. 

the  square  domain  D.  The  same  arguments  allow  us  to  write 


(0V20  +  0V20)  = 


</>VV  =  V  •  (0V(V2</>)  -  (V0)(V20))  +  |V20|2. 

Collecting  all  the  above  results,  one  can  write 

|  X  W  =  2I JD  -  Id2)  +  *uJD  |V0|2  -  <r2  /  IVVI2],  (21) 

and  also 

jt  \<t>?  =  2[ j  |</>|2(a  +  w2  -  |0|2)  -  J  | eV20  +  cu0|2],  (22) 

as  it  can  be  readily  seen  after  taking  into  account  the  following  relation 

|  eX72(j)  +  co(j)\2  =  e2\V2(j)\2  +  u2\(t)\ 2  +  £ou$V2(j)  +  0V20) 
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and  the  result  above  for  the  integral  of  its  last  term. 

Finally,  for  u  <  0  all  negative  terms  in  the  right  hand  side  of  eq.  (21)  can  be 
removed  to  give 

iLM2<-2aL'^ 

that,  according  to  Gronwall’s  Lemma,  implies  that 


<  Ce 


2  at 


0  as  t  — >  oo, 


id 


and  thus  we  can  conclude  that 


/  \(p\J  — >  0  as  t  — >  oo, 

J  D 

i.e.,  all  solutions  decay  to  zero  for  a  <  0  and  u  <  0.  If  now  u  >  0,  then  we  can  use 
Eq.  (22)  and  follow  the  same  argumentation  as  above  to  obtain  that,  for  a  +  co2  <  0, 
all  solutions  also  decay  to  zero. 

It  is  interesting  to  notice  that  the  global  stability  of  the  zero  solution  allow  us  to 
conclude  that  the  CSH  cannot  exhibit  localized  solutions  with  a  nonzero  central  core 
and  tail  that  decays  to  zero;  in  order  to  have  such  solutions  the  zero  solution  has 
to  be  stable  but  this  also  implies,  because  of  the  global  stability,  that  the  pulse-like 
solution  has  to  decay  to  zero. 


4.1  TW  solutions 

The  simplest  nonzero  solutions  of  the  CGL  eq.  (18)  are  the  spatially  uniform  trav¬ 
elling  waves  (TW), 

0tw  =  eikTw-x_ifcTwt)  (23) 

with  /cTW  =  |kTW|  ~  1,  which  exist  only  for  a  >  0  and  are  approximate  solutions  of 
the  CSH  up  to  0(e)  corrections. 

The  TW  stability  characteristics  can  be  easily  obtained  looking  for  solutions  of 
the  form 

0  =  0tw(1  +  0  with  |£|  <C  1 

in  the  CSH  and  neglecting  all  but  the  linear  terms.  The  following  equation  for  the 
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Figure  6:  Growth  rate  of  the  perturbations  of  the  TW  for  ui  >  0  and  lo  <  0,  with  the  two  wavenum¬ 
ber  regions  shown:  k  ~  1  (dispersive  perturbations)  and  k  ~  1  j yfe  1  (diffusive  perturbations). 

evolution  of  the  perturbation  £  is  obtained 

6  =  i(V2£  +  i2kTW  •  V£)  -  a(£  +  £)  -  2au[V2£  +  ...]-  0[V4£  +  ...], 

where,  as  we  will  see  below,  only  the  main  contribution  of  the  small  diffusive  terms 
needs  to  be  retained.  If  we  now  expand  £  in  Fourier  series 

«  =  £&Me‘k'x. 

only  the  mode  pairs  with  wavevectors  ±k  are  coupled  and  their  evolution  equations 
read 

=  (—a  —  i k2  —  i2kTW  •  k  +  2 etok2  —  £20)£k  — 

— =  (—a  +  i k2  —  i2kTw  •  k  +  2eujk 2  —  £2/c4)£_k  —  cc£k- 

where  k  —  |k|.  The  solutions  of  this  linear  system  are  proportional  to  ext  with  A 
given  by 

A  =  —a  —  i2kTw  •  k  +  2eu :k2  —  £2kA  ±  V a2  —  k4. 

Again,  we  can  distinguish  between  perturbations  exhibiting  only  dispersive  scales, 
i.e.,  with  wavevector  |k|  ~  1,  which  have  a  growth  rate  that  can  be  approximated 
as 

Re(A)  =  -a±  Va2  -  kA  +  . . .  (24) 
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and  are  always  stable  because  a  is  positive  for  all  TW  (see  Fig.  6);  this  is  again  the 
result  that  is  obtained  if  the  linear  stability  analysis  is  performed  on  the  CGL.  And 
perturbations  with  high  diffusive  wavevectors,  k  ~  K/V?  1  ,  whose  dynamics  is 
not  contained  in  the  CGL,  and  their  growth  rate,  up  to  0(y/e)  corrections,  is  given 
by 

Re(A)  =  -a  +  2uK2  -  K 4  +  . . . ,  (25) 

which  is  also  plotted  in  Fig.  6.  These  diffusive,  short  wave  perturbations  become 
unstable  outside  the  region  defined  by  u  >  0  and  a  >  oj2  (marked  with  light  shading 
in  Fig.  5)  with  critical  wavenumber  k  =  \/uj \fe  3>  1.  Note  that  the  growth  rate  of 
the  diffusive  perturbations  is  of  order  one  and  thus  this  not  a  higher  order  correction 
that  evolves  on  a  much  slower  time  scale;  this  instability  develops  in  the  time  scale 
t  ~  1  of  the  dispersive  CGL  dynamics. 

Finally,  it  is  also  interesting  to  mention  that  the  region  u  >  0  and  —  u2  <  a  <  u2 
in  Fig.  5  appears  to  be  the  best  suited  to  look  for  complex  dynamics  because  in 
there  both  the  zero  solution  and  the  uniform  TW  are  unstable. 

4.2  Spirals 

The  rotating  spiral  wave  solutions  of  the  CGL  are  of  the  form 

0s  =  yWp(r)  eH"S t+me+i,{r))^ 

where  (r,  6)  are  polar  coordinates,  m  e  N  represents  the  number  of  arms  of  the  spiral 
and  cns  its  rotating  angular  frequency.  After  inserting  the  expression  above  into  the 
CGL  eq.  (18),  the  following  equations  are  obtained  for  the  functions  p(r)  and  0(r) 

|  'TTl2 

Prr  +  ~Pr  ~  p{^r  +  ~^)  =  ^SP,  (26) 

4>rr  +  “tV  +  2— =  1  -  p2,  (27) 

r  p 

together  with  the  following  boundary  conditions  that  prevent  the  formation  of  sin¬ 
gularities  at  r  =  0 

0r( 0)  =  0,  p(r)  =  p0rm  for  r  — >  0,  and  p(r)  — >  1  for  r  — >  oo,  (28) 
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see  e.g.  [24],  The  spiral  is  zero  at  the  origin  and,  as  r  — >  oo,  its  amplitude  approaches 
that  of  the  uniform  TW  y/a  and  its  radial  wavenumber  verihes 

A  ->•  Aooi  With  US  =  -Iproo, 


which  is  precisely  the  frequency-wavenumber  relation  of  the  TW.  These  spirals  of 
the  CGL  are  again  solutions  of  the  CSH  up  to  0(e)  error. 

The  numerical  integration  of  the  semi-infinite,  singular  boundary  value  problem 
(26)- (28)  has  been  performed  following  the  steps  below. 


1.  In  order  to  avoid  the  singularity  at  r  =  0,  the  solution  at  r  —  r0  <C  1  is 
approximated  by  the  expansions  for  r  — >  0 


P(r) 


=  p0rm  + 


AWS  ^,m+ 1  _i_ 

4  (m  +  1) 


A  (r) 


r 

2  (m  +  1) 


+  .... 
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2.  The  third  order  system  (26)-(27)  for  (pr,p,  A)  is  now  integrated  from  r  =  r0 
to  r  =  Too  S>  1  and,  at  r  =  r^,  the  solution  is  forced  to  match  the  expansion 
for  r  — >  oo 

=  _  ^Voo  _  _  4m2  +  29^  _  80m2  +  973^  +  32 

P[r)  2  r  8  r2  16^roor3  128r4 

3.  For  every  fixed  value  of  m,  the  matching  conditions  for  p  and  pr  at  r  =  Too 
give  of  a  two  equation  system  for  us  and  p0  that  is  solved  applying  Newton’s 
method  (with  the  Jacobian  approximated  using  second  order  finite  differences). 
Note  that  there  is  no  need  to  match  since  it  is  completely  defined,  once  p  is 
known,  by  the  expression 

i  r 

Aij)  =  /  sp2(s)(l  -  p2(s))  ds, 

rp\r)  J o 

which  is  readily  obtained  from  eq.  (27)  after  multiplying  by  rp 2  and  integrating 
from  0  to  r. 
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Ill  —  1 


Figure  7:  From  top  to  bottom:  one,  two  and  three  armed  spiral  solutions  corresponding  to  in  = 
1,2  and  3.  The  left  column  shows  the  amplitude  p  and  radial  wavenumber  ij)r  profiles.  Right 
column  shows  the  grayscale  map  of  the  real  part  of  <j> g. 

The  resulting  spirals  for  m  =  1,2  and  3  are  shown  in  in  Fig.  7,  and  the  corre¬ 
sponding  values  of  po  and  the  rotating  frequency  cus  are  given  in  the  table  below. 


m 

us 

Po 

i 

-0.426468 

0.365958 

2 

-0.700454 

0.0957711 

3 

-0.909445 

0.0192913 

20 


The  computations  have  been  done  for  tq  =  1CT3  and  =  20  in  order  to  ensure 
that  the  errors  in  the  expansions  for  r  — >  0  and  r  — >  oo  above,  which  are  ~ 
and  ~  I/t^q  respectively,  are  both  below  10~6  (we  have  also  checked  that  the  digits 
presented  in  the  table  do  not  change  when  r0  is  reduced  to  10~4  and  r <*,  is  increased 
to  30). 
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5  Numerical  results 


The  purpose  of  the  numerical  results  presented  in  this  section  is  twofold:  (i)  We 
will  first  perform  some  numerical  simulations  of  the  order  parameter  equation  (the 
CSH  eq.  (17))  to  check  the  stability  predictions  of  the  previous  section  and  to  show, 
as  it  was  advanced  in  Section  3,  the  presence  of  two  essentially  different  scales  in 
the  extended  system  limit  £  — »  0,  namely,  dispersive,  <5disp  ~  1,  and  diffusive  scales, 
<5diff  ~  And  (ii)  we  will  then  numerically  integrate  both,  the  complete  MB  eqs. 
(l)-(3)  and  the  CSH  eq.  (17),  and  compare  their  respective  solutions  to  confirm  the 
validity  of  the  CSH  as  the  appropriate  order  parameter  equation  for  the  description 
of  the  weakly  nonlinear  dynamics  of  the  MB. 

5.1  Numerical  simulations  of  the  CSH 

The  CSH  equation  (17)  is  numerically  integrated  in  a  square  domain  [0, 1]  x  [0.1] 
with  periodic  boundary  conditions  for  four  representative  sets  of  parameter  values. 


a 


Figure  8:  Stability  diagram  for  the  zero  solution  and  the  TW.  Light  shading:  TW  stable.  Dark 
shading:  zero  solution  stable.  Numbered  dots  correspond  to  the  parameter  values  used  in  the 
numerical  integrations. 


CASE  1:  a  =  —0.75  and  ca  =  0.5.  For  these  parameter  values,  the  results 
obtained  in  Section  4  indicate  that  the  zero  solution  is  globally  stable,  see  Fig.  8. 
The  initial  condition  used  in  all  simulations  is  a  constant  amplitude  pure  mode  with 
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wavevector  kxw  =  (1,2),  represented  in  Fig.  9,  and  with  a  superimposed  random 
perturbation  of  amplitude  10~3.  The  CSH  has  been  integrated  for  three  different 
values  of  £  and  the  resulting  time  evolution  of  the  norms  of  the  solution 


1101 


\4>\dxdy  and  ||V0|| 


1 1  V0|  |  dx  dy 


are  shown  in  Fig.  10.  As  it  was  expected,  the  dynamics  of  the  system  in  this  case 
is  extremely  simple,  and  all  three  solutions  decay  to  zero.  It  is  interesting  to  notice 
that  the  decay  rate  for  £  -C  1  appears  to  be  independent  of  £  (see  the  slopes  of  the 
two  solid  lines  in  the  lower  plot  in  Fig.  10)  but  different  from  that  for  e  =  0  (dashed 
line  in  the  lower  plot  in  Fig.  10),  which  decays  faster.  This  is  also  in  agreement  with 
the  results  from  Section  4:  the  dispersion  relation  of  the  zero  solution  for  e  1  is 
given  by  eq.  (20),  whose  maximum  value  a  —  to2  =  —0.5  (associated  to  the  slowest 
decaying  modes)  is  independent  of  e  and  different  form  that  corresponding  to  e  =  0 
(when  the  CSH  equation  simplifies  to  the  the  CGL  equation  (18))  that,  according 
to  eq.  (19),  is  constant  and  equal  to  a  =  —0.75,  and  thus  produces  a  faster  decay 
to  zero. 


Re(</>)  Im  ((/)) 


Figure  9:  Colored  contour  maps  of  the  real  (left)  and  imaginary  (right)  part  of  a  pure  mode  field 
of  the  form  (j>  =  y/[ofelk'x  with  a  =  —0.75  and  kxw  =  (1,2).  Lighter  (darker)  colors  correspond 
to  higher  (lower)  values. 
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1 


Figure  10:  Time  evolution  of  the  norm  of  the  solution  of  the  CSH  eq.  (17)  (above)  and  its  gradient 
(below).  The  initial  condition  is  the  pure  mode  in  Fig  9  with  a  random  perturbation  of  size  10-3, 
and  the  parameters  correspond  to  the  CASE  1  in  Fig.  8  with  e  =  10-3  (thick  line),  e  =  10~3/4 
(thin  line)  and  e  =  0  (CGL  eq.  (18),  dashed  line). 

CASE  2:  a  =  0.75  and  u  =  0.5.  We  are  now  in  the  region  of  the  parameter 
space  where  the  uniform  TW  described  by  expression  (23)  in  Section  4.1  are  all 
stable  (see  Fig.  8).  Recall  that  these  are  TW  with  wavenumber  k  ~  1  (without 
diffusive  scales)  and  therefore,  in  the  £<1  limit,  they  are  approximately  given  by 
the  CGL  equation  (18).  The  starting  point  of  the  numerical  integrations  shown  in 
Fig.  11  is  now  a  uniform  TW  (i.e.  that  with  kxw  —  (0,0)  and  amplitude  ^/a)  with 
a  10~3  random  perturbation,  and,  as  expected,  the  system  relaxes  to  the  uniform 
TW:  1 10||  —>  y/a  and  ||V0||  —>  0.  As  it  can  be  appreciated  from  the  lower  plot  of 
Fig.  11,  the  resulting  decay  rate  for  e  -C  1  is  again  independent  of  £  and  is  given  by 
the  maximum  of  the  growth  rate  (25),  —a  +  lu2  =  —0.5,  and  it  again  differs  from 
that  for  e  =  0  that,  according  to  expression  (24),  is  equal  to  —a  =  —0.75. 
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Figure  11:  Time  evolution  of  the  norm  of  the  solution  of  the  CSH  eq.  (17)  (above)  and  its  gradient 
(below).  The  initial  condition  is  a  spatially  uniform  steady  solution  with  a  random  perturbation 
of  size  10-3,  and  the  parameters  correspond  to  the  CASE  2  in  Fig.  8  with  e  =  10-3  (thick  line), 
e  =  10_3/4  (thin  line)  and  e  =  0  (CGL  eq.  (18)  dashed  line). 


Figure  12:  Time  evolution  of  the  norm  of  the  solution  of  the  CSH  eq.  (17)  (above)  and  its  gradient 
(below).  The  initial  condition  is  a  spatially  uniform  steady  solution  with  a  random  perturbation 
of  size  10-3,  and  the  parameters  correspond  to  the  CASE  3  in  Fig.  8  with  e  =  10~3  (thick  line), 
e  =  10_3/4  (thin  line)  and  e  =  0  (CGL  eq.  (18),  dashed  line). 
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CASE  3:  a  =  —0.5  and  u  —  2.  The  initial  condition  is  the  same  perturbed 
uniform  state  used  in  CASE  2,  but  now,  as  it  is  clear  from  Fig.  12,  the  behaviour 
of  the  system  for  e  <C  1  is  completely  different  from  that  for  e  =  0.  The  resulting 
temporal  evolution  for  e  =  0  is  just  a  monotonous  decay  to  the  zero  solution  (dashed 
line  in  Fig.  12)  that,  according  to  eq.  (19),  is  always  stable  in  the  CGL  dynamics 
when  a  is  negative.  On  the  other  hand,  the  zero  solution  is  unstable  for  the  CSH 
with  e  -C  1  (see  Fig.  8),  and  what  happens  now  is  that  the  system  evolves  to  a 
nonzero  final  state  that  exhibits  small  diffusive  scales  (solid  lines  in  Fig.  12).  These 
small  diffusive  scales  grow  exponentially  with  a  finite  non-zero  growth  rate  as  e  — >  0, 
see  the  lower  plot  of  Fig.  12.  Note  that,  despite  of  the  small  coefficients  of  diffusion 
and  double  diffusion  in  the  CSH  eq.  (17),  the  onset  of  these  small  diffusive  scales  is 
not  a  higher  order,  longer  time  scale  effect;  they  evolve  in  the  t  ~  1  time  scale  of  the 
CSH,  as  it  is  clear  from  the  slope  of  the  solid  lines  in  the  lower  plot  of  Fig.  12.  The 
resulting  values  of  ||V0||  at  t  =  20  for  e  =  10”3  and  £  =  10-3/4  are,  respectively, 
82.187  and  164.85,  whose  ratio  is  equal  to  2.006  . . .  ,  in  agreement  with  the  fact  that 
the  small  diffusive  scales  dominate  the  final  state,  forcing  the  norm  of  the  gradient 
to  behave,  in  first  approximation,  as  ||V0||  ~  l/y/£.  Moreover,  the  solutions  at 
t  =  20  for  £  =  10~3  and  £  =  10~3/4  (shown  in  Figs.  13  and  14)  resemble  a  TW  with 
a  number  of  wavelengths  that  is  roughly  7  and  14,  indicating  that  this  final  state 
corresponds  to  a  TW  but  with  a  diffusive  wavenumber  k  1/VS- 


Re(0)  lm(0) 


Figure  13:  Colored  contour  maps  of  the  real  (left)  and  imaginary  (right)  part  of  the  solution  of  the 
CSH  eq.  (17)  at  t  =  20  for  a  =  —0.5,  u  =  2  and  £  =  10-3. 
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Figure  14:  Colored  contour  maps  of  the  real  (left)  and  imaginary  (right)  part  of  the  solution  of  the 
CSH  eq.  (17)  at  t  =  20  for  a  =  -0.5,  uj  =  2  and  e  =  10~3/4. 


CASE  4:  a  =  0.5  and  u  —  2.  The  initial  condition  is  now  that  of  the  CASE 
1  simulations:  a  TW  with  wavevector  kxw  —  (1)2)  shown  in  Fig.  9  with  a  10~3 
perturbation.  This  TW  is  stable  in  the  CGL  dynamics;  this  is  clearly  seen  in  the 
results  for  e  =  0  plotted  in  Fig.  15,  where  the  norm  of  the  solution  remains  almost 
constant  and  equal  to  its  initial  value  y/a,  and  the  perturbations  (lower  plot)  decay 
exponentially.  As  it  happened  in  the  previous  case,  this  TW  is  unstable  for  the  CSH 
with  e<l,  and  the  solution  develops  small  diffusive  scales  that  grow  exponentially 
(see  the  solid  lines  in  Fig.  15).  The  ratio  of  the  final  values  of  ||V0||  for  £  =  10~3 
and  e  =  lCC3/4  is  approximately  176.2/87.6  =  2.01,  and  there  are  twice  as  much 
wavelengths  in  the  solution  in  Fig.  17  for  e  =  10~3/4  than  in  that  shown  in  Fig.  16 
for  e  =  10~3,  indicating  again  that  the  final  state  is  a  TW  with  diffusive  wavenumber 
k  ~  1  j \fe.  Notice  that  this  TW  is  very  similar  to  that  obtained  in  the  previous  case; 
this  is  due  to  the  fact  that  the  most  unstable  wavenumber  is  given  by  k  ~  y/ui/y/e 
(see  eq.(25))  that  depends  only  on  u  and  £,  which  take  the  same  values  in  both  cases 
(the  amplitude  of  both  TW  is  different  but  this  cannot  be  appreciated  in  the  color 
maps) . 
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Figure  15:  Time  evolution  of  the  norm  of  the  solution  of  the  CSH  eq.  (17)  (above)  and  its  gradient 
(below).  The  initial  condition  is  a  spatially  uniform  steady  solution  with  a  random  perturbation 
of  size  10-3,  and  the  parameters  correspond  to  the  CASE  4  in  Fig.  8  with  e  =  10-3  (thick  line), 
e  =  10_3/4  (thin  line)  and  e  =  0  (CGL  eq.  (18),  dashed  line). 


Re(0)  Im  (0) 


Figure  16:  Colored  contour  maps  of  the  real  (left)  and  imaginary  (right)  part  of  the  solution  of  the 
CSH  eq.  (17)  at  t  =  20  for  a  =  0.5,  u  =  2  and  e  =  10-3. 
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Figure  17:  Colored  contour  maps  of  the  real  (left)  and  imaginary  (right)  part  of  the  solution  of  the 
CSH  eq.  (17)  at  t  =  20  for  a  =  0.5,  w  =  2  and  e  =  lCT3/4. 


The  numerical  simulations  presented  above  clearly  show  that  there  are  two  types 
of  solutions  of  the  CSH  with  a  completely  different  behaviour  in  the  physically 
relevant  limit  of  large  systems  £  — >  0:  those  that  only  exhibit  scales  that  are  of 
the  order  of  the  total  size  of  the  domain,  5disp  ~  1,  and  that  can  be  accurately 
approximated  by  the  simpler  CGL  that  is  obtained  by  setting  e  =  0  in  the  CSH 
(CASES  1  and  2),  and  those  that  develop  small  diffusive  scales,  <5difr  ~  y/s,  whose 
dynamics  can  only  be  correctly  described  using  the  CSH  with  £<1  (CASES  3  and 

4)- 

The  details  of  the  applied  numerical  scheme  are  given  below  in  section  5.4.  The 
typical  number  of  Fourier  modes  used  in  the  simulations  above  are  64  x  64  and 
128  x  128,  with  a  time  step  in  the  range  dt  =  0.001 . . .  0.0001. 

5.2  MB  and  CSH  solution  comparison 

In  this  section  we  numerically  check  the  accuracy  of  the  CSH  eq.  (17)  as  the  order 
parameter  equation  for  the  description  of  the  weakly  nonlinear  dynamics  of  the  MB 
eqs.  (l)-(3)  with  small  detuning. 

We  numerically  integrate  the  MB  in  a  large  square  domain  [0,  L\  x  [0,  L\  with 
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periodic  boundary  conditions  and  L  >  1.  The  pump  and  the  detuning  are  taken 
appropriately  small 

|r  —  1|  ~  —  -C  1  and  |Q|  ~  —  <C  1, 

1J  ±J 

and  the  initial  condition,  according  to  the  scaling  in  (16),  takes  the  form 

1 
1 
0 

where  b  and  cr  are  fixed,  order  one  parameters  (class  C  laser),  and  the  complex 
function  0O  is  periodic  in  [0, 1]  x  [0, 1]  and  such  that  |0O|  ~  1;  and  the  final  integration 
time  fMB  is  chosen  such  that 

£mb  ~  L2  3>  1. 

We  also  integrate  the  CSH  in  the  domain  [0, 1]  x  [0, 1]  with  periodic  boundary 
conditions,  4>o(x,y )  as  initial  condition  and  the  following  parameters  (see  (16)) 


(29) 

(30) 

(31) 


and  up  to  a  final  time 


1 

^CSH  —  7 - ,  iUJmB' 

(a  +  1)L2 


(32) 


In  order  to  quantitatively  measure  the  accuracy  of  the  order  parameter  descrip¬ 
tion,  we  first  perform  several  simulations  of  the  MB  with  increasing  domain  size  and 
final  time 
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Figure  18:  Colored  contour  maps  of  the  real  (left)  and  imaginary  (right)  part  of  the  initial  condition 
used  for  the  CSH  and  MB  comparison. 


and  the  corresponding  ones  of  the  CSH 


CSHi 

csh2 

csh3 

£l 

£2 

£3 

keeping  fixed  b ,  a ,  a,  uj  and  the  final  CSH  integration  time,  tcsu,  and  obtaining  the 
rest  of  the  parameters  (rt  —  1),  £mb;  and  e*  from  the  above  eqs.  (29)-(32)  for 
every  given  value  of  Lt.  And  then,  we  evaluate  the  difference  of  both  solutions 


di(t)  = 


E 

P 

N 


a 


-I  MB, 


1 

1 

0 


(j>csHi(x / Li,  y/ Li,  tcsKj)  ^ 


and  check  the  asymptotic  behavior  predicted  in  section  3, 

di/£i  — >  0  as  £i  — »■  0. 


(34) 


We  repeat  the  above  procedure  for  two  sets  of  parameter  values: 

CASE  1:  a  —  1,  b  —  1,  a  —  0.75  and  u  =  0.5.  The  initial  condition  (plotted  in 
Fig.  18)  is  given  by  the  expression 

<j>o{x,y)  =  0.5  +  3i  +  cos(x  +  cosy)esinx+i3cos2y, 
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which  presents  several  peaks  and  valleys  and  has  been  selected  in  order  to  ensure  that 
the  simulations  are  started  at  a  generic  state,  away  from  the  simple  TW  solutions. 
Both  systems  have  been  integrated  up  to  final  time  fcsH  =  10  and  for  four  different 
values  of  L ,  with  parameters 


And  the  resulting  MB-CSH  differences  (33)  are  shown  in  Fig.  19.  The  time  evolution 
of  the  MB  presents  a  short  initial  stage,  that  takes  place  for  tMB  ~  1  (  tcsH  ~ 
1  / L2  <C  1),  in  which  the  system  quickly  damps  out  all  modes  with  high  wavenumber. 
This  produces  a  sharp  growth  and  decrease  of  the  MB-CSH  difference  that  can  be 
appreciated  in  the  lower  left  plot  of  Fig.  19.  After  this  initial  phase,  the  system 
enters  the  slow  dynamics  regime  and  its  evolution  takes  place  in  the  scale  fcsH  ~  1 
(^csh  ~  L2  S>  1),  see  top  plot  of  Fig.  19.  The  maximum  values  of  each  difference, 
dmax)  divided  by  £  are  plotted  against  e  in  the  lower  right  plot  of  Fig.  19;  it  is  clear 
from  this  plot  that  dm^/s  goes  to  zero  as  £  is  decreased,  confirming  the  asymptotic 
behavior  predicted  in  (34). 

CASE  2:  a  —  1,  b  —  1,  a  —  0.5  and  u  =  2.  The  initial  condition  is  that  of 
the  previous  case  and,  again,  both  systems  have  been  integrated  up  to  final  time 
tcsH  =  10,  but  now  for  five  different  values  of  L,  with  parameters 


and  the  results  are  now  presented  in  Fig.  20.  Again  the  system  exhibits  a  fast 
adaptation  stage  (lower  left  plot  of  Fig.  20)  followed  by  a  slow  evolution  in  the 


32 


tcsH  ~  1  time  scale  (top  plot  of  Fig.  20).  But  now,  we  are  in  a  point  of  the 
parameter  space  where  the  final  state  exhibits  diffusive  scales  hdifr  ~  VL  in  the  MB 
variables.  And  thus,  it  is  now  necessary  to  go  to  higher  L  values  (L  >  50,  £  <  .01)  in 
order  to  have  a  clear  separation  of  scales  between  <5disp  ~  L  and  <5diff  ~  VL  (hdisp  ~  1 
and  hdiff  ~  y/s  hr  the  CSH  scaling)  and  start  to  see  the  asymptotic  decay  predicted 
in  (34),  see  the  dots  in  the  bottom  right  plot  of  Fig.  20. 

The  numerical  method  used  to  integrate  the  MB  eqs.  (l)-(3)  is  briefly  described 
in  the  following  section.  The  computations  above  have  been  carried  out  with  64  x  64 
Fourier  modes  and  time  integration  step  dt  =  0.0025  and  dt  =  0.00125.  Note  that 
such  small  time  steps  are  necessary  in  order  to  be  sure  that  the  numerical  integration 
errors  of  the  MB  equations  are  kept  small  and  do  not  contaminate  the  differences 
that  we  trying  to  compute.  The  combination  of  small  time  steps  (dt  =  0.0025) 
and  very  long  integration  times  (£mb  =  50000)  makes  the  computations  presented 
above  extremely  CPU  costly;  the  computation  of  the  time  evolution  of  one  MB-CSH 
difference  may  take  up  to  5  days  in  a  3GHz  Pentium  IV  computer  with  1Gb  RAM. 
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Figure  19:  Top  plot:  Time  evolution  of  the  difference  of  the  solution  of  the  CSH  and  MB  equations 
for  a  =  1,  b  =  1,  a  =  0.75,  u>  =  0.5  and  L  =  20,30,40  and  50.  Bottom  plot:  short  time  detail  of 
evolution  of  difference  (left)  and  (right)  the  resulting  maximum  values  of  the  difference  versus  e 
(the  dashed  line  corresponds  to  a  decay  ~  e  ). 
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Figure  20:  Top  plot:  Time  evolution  of  the  difference  of  the  solution  of  the  CSH  and  MB  equations 
for  a  =  1,  b  =  1,  a  =  0.5,  u>  =  2  and  L  =  20,30,40,50  and  60.  Bottom  plot:  short  time  detail  of 
evolution  of  difference  (left)  and  (right)  the  resulting  maximum  values  of  the  difference  versus  e 
(the  dashed  line  corresponds  to  a  decay  ~  e  ). 


35 


5.3  Numerical  integration  of  the  MB  equations 

The  Maxwell-Bloch  equations  (l)-(3)  with  periodic  boundary  conditions  in  a  square 
domain  of  size  L , 

(E,  P,  N)  (x  +  L,  y,t)  =  (E,  P,  N)  (x,  y,t) ,  (E,  P,  N)  (x,  y  +  L,t)  =  (E,  P,  N)  ( x ,  y,  t ) , 

are  numerically  integrated  using  a  Fourier  series  representation  in  x  and  y,  and  a 
4th  order  Runge-Kutta  temporal  scheme. 

In  order  to  do  this,  we  expand  the  solutions  in  a  double  Fourier  series  of  the  form 

( E,P,N)(x,y,t )  =  ^(ekj(t),Pkj(t),nkj(t))e^2v/L)k'x, 

k,j 

with  k  =  (k,j)  and  x  =  (x,y),  and 


^k.j  Tl—k,—ji 

because  the  population  inversion  N(x,  y,  t )  is  a  real  valued  field.  The  resulting 
equations  for  the  Fourier  components  can  be  written  as 

<^L  =  -(<7  +  i(27r /L)2(k2  +  j2))ekJ  +  apkJ: 

=  -(1  +  i El)pk,j  +  rekJ  -  [NE]ktj, 

^  =  -tokJ  +  l[EP  +  EP]kJ. 

We  have  to  use  a  sufficiently  large  number  of  Fourier  modes  in  order  to  ensure  that 
we  are  accurately  solving  the  shortest  spatial  scales  exhibited  by  the  system,  which 
are  A  ~  \/~L.  This  implies  that  the  number  of  Fourier  modes,  should  be  such 
that 

Nf  rsj  m\\[L , 

where  rri\  represents  the  number  of  modes  used  to  describe  the  shortest  scales,  and 
m\  ~  5  —  10  typically.  Therefore,  the  coefficient  that  comes  from  dispersion  in  the 
system  above  can  be  quite  large  and  this  can  give  rise  to  numerical  instabilities  in 
the  usual  explicit  time  integration  schemes  (see  [25]). 
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This  problem  can  be  avoided  if  we  integrate  explicitly  the  linear  diagonal  terms 
to  obtain 


e~Clk’jtUik,j 

e~Clkjt flkj(uik> J>,  U2k> j>,  U3kf  J>) 

e~C2k’jtU2k,j 

= 

e-c2k,jt  f2k  j(u\k, J,,  U2k’j',  U3k‘ J/) 

_  e~C3k’jtu3kJ  _ 

_  e_C3fc^*/3jfcj(^ife,j',  u2k>tfi  u3k>tf) 

where 


^ k,j  (iti,  1l2i  U3)kj  {eiPi  n)k,j, 

ck,j  =  (ci,c2,c3)fej  =  (-(a  +  i(27r/L)2(/c2+j2)),-(l  +  ifi),-6), 
fkj  =  (A,  /2,  fz)k,j  =  (vu2k,j,  rulkJ  -  [NE\kJ,  -[ EP  +  EP]kJ), 

And  then,  we  use  a  standard  explicit  fourth  order  Runge-Kutta  method  [25],  which 
yields  the  following  time  marching  scheme  for  a  single  time  step  of  size  At 


uktj{t  +  A  t)=  ec^Atukj(t)+  f  [e^MKld  +  2ec^At/2K'lj+  2ec^At/2K3kd+KlJ], 

where  the  vectors  Kkj-  Kk-  and  Kk-  are  given  by 

^k,j  —  fkj{uk'j'(t ))) 

Klj  =  fk,j  (eCk'j'At/2Uk'f  (t)  + 

Kl,j  =  fk,j  (eCk'j,At/2uk'j'  (t)  +  YKh')’ 

K£j  =  fkj(eck'i'Atuklf(t)  +  A  teW^Klj'), 

and  the  products  of  the  vectors  by  the  exponentials  of  the  form  eCk-jAt  have  to  be 
performed  component  by  component. 

This  integration  procedure  always  handles  the  Fourier  representation  of  the  so¬ 
lution  and  it  is  transformed  to  physical  space  only  during  the  computation  of  the 
nonlinear  terms,  which  is  performed  in  the  following  steps 

ukJ  tQ  PhyS,>  (E,  P,  N)  -  (aP,  rE  -  NE,  l-[EP  +  EP])  --  ---  fkj, 
where  the  products  of  the  amplitudes  are  computed  in  physical  space  and  the  so- 
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called  2/3  rule  is  used  to  remove  the  aliasing  terms  (see  e.g.  [26]). 

This  numerical  procedure  has  been  implemented  in  a  FORTRAN  code  and  the 
FFTW  subroutines  [27]  are  used  to  perform  the  Fourier  transforms. 


5.4  Numerical  integration  of  the  CSH  equation 

The  numerical  integration  method  for  the  CSH  equation  (17)  in  the  square  domain 
[0.1]  x  [0, 1]  with  periodic  boundary  conditions, 

4>{x  +  1  ,y,  t )  =  0(x, y,  t),  y  +  l,t)  =  y,  t), 


is  completely  similar  to  that  described  in  the  previous  section  for  the  MB.  First  the 
solution  is  expanded  in  Fourier  series 

<f>(x,y,t)  =  5^0fcj(t)el(2,r)k'x, 

k,j 

with  k  =  (k,j)  and  x  =  (x,y),  and  then  the  resulting  system  of  ODE’s  for  the 
Fourier  modes 

=  (a  -  (2n)2(k2  +  j2){2euj  +  i)  -  £2(2n)4(k2  +  j2)2)(pkJ  -  [</#[ 2}kj 

is  integrated  with  the  same  temporal  marching  scheme  used  for  the  MB.  The  linear 
terms  are  integrated  explicitly  with 

Cfcj  =  a  —  (27r)2(/c2  +  j2)(2euj  +  i)  —  £2(2n)4(k2  +  j2)2,  and 
fkj  =  [0101 2]k,j, 


and  the  nonlinear  terms  are  again  computed  in  physical  space 


to  Phys. 

rk,j  * 


to  Fou. 


/fc. 


,3i 


using  the  2/3  rule  for  the  aliasing  terms  (see  e.g.  [26]).  In  this  case,  the  smallest 
spatial  scale  that  we  have  to  represent  are  the  diffusive  scales  <5diff  ~  -C  1 
and  thus  the  required  number  of  Fourier  modes  goes  as  N-p  ~  rris / \fe  3>  1  with 
m§  ~  5  —  10  typically. 
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This  numerical  procedure  has  been  implemented  in  a  FORTRAN  code  and  the 
FFTW  subroutines  [27]  are  used  to  perform  the  Fourier  transforms.  This  code  is 
also  used  to  integrate  the  CGL  eq.  (18)  just  by  setting  £  =  0. 
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